library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6     ✔ purrr   0.3.4
✔ tibble  3.1.7     ✔ dplyr   1.0.9
✔ tidyr   1.2.0     ✔ stringr 1.4.0
✔ readr   2.1.2     ✔ forcats 0.5.1
── Conflicts ───────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
library(plotly)
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

Consider a somewhat nasty-looking curve in two dimensions. We define \(z\) to be a function of \(x\) and \(y\) that looks like this: \[ z = f(x, y) = \biggl[1 + e ^{-|x + 2y|} \cos\biggl(0.2 \sqrt{x^2 + y^2}\biggr) \biggr] \times I\{-10 < x,y < 10 \} \] In other words, it is a function that returns 0 everywhere outside of a square of width 20, centered on the origin. But within that square it returns positive values.

Here is some code that defines it:

ugly <- function(x, y) {
  ifelse(x < -10 | y > 10 | x > 10 | y < -10, 
         0,
         1 + exp(-abs(x + 2 * y) / 10) * cos(0.2 * sqrt(x^2 * y^2))
  )
}

And, we can use that to visualize the surface by evaluating a lot of points and plotting it with the plotly package. Go ahead and rotate and manipulate the thing.

vals <- seq(-12, 12, by = 0.05)
pts <- expand_grid(
  x = vals,
  y = vals
) %>%
  mutate(z = ugly(x, y))

zmat <- matrix(pts$z, ncol = length(vals), nrow = length(vals))

plot_ly(x = vals, y = vals, z = ~zmat) %>%
  add_surface()

The volume under that curve is given by the definite integral: \[ \int_{-10}^{10} \int_{-10}^{10} f(x, y) dxdy = \int_{-10}^{10} \int_{-10}^{10} \biggl[1 + e ^{-|x + 2y|} \cos\biggl(0.2 \sqrt{x^2 + y^2}\biggr) \biggr] dxdy \] That looks a little hairy. It may have an analytical solution, or it could be evaluated numerically with high accuracy. But, if you rewrite that integral as an expection over a two-dimensional uniform variable on the square from -10 to 10 in \(x\) and also from -10 to 10 in \(y\), then you get: \[ 20 \times 20 \times \int_{-10}^{10} \int_{-10}^{10} f(x, y) \frac{1}{20}\frac{1}{20}dxdy = 400 \times \mathbb{E}[f(x, y)] \] where the expectation is taken over a 2-D uniform variable from -10 to 10 in both x and y.

Given that as a preamble, your mission it to approximate that integral using Monte Carlo sampling.

Good luck! The answer is \(\approx 423.7\).

Here is how that can be done:

set.seed(123)
xs <- runif(1e6, -10, 10)
ys <- runif(1e6, -10, 10)

answer <-  400 * mean(ugly(xs, ys))

answer
[1] 423.7917
LS0tCnRpdGxlOiAiQm9udXMgZXhjZXJjaXNlOiB1c2UgTW9udGUgQ2FybG8gdG8gYXBwcm94aW1hdGUgYSBkZWZpbml0ZSBpbnRlZ3JhbCAtLSBLRVkiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCgpgYGB7ciwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkocGxvdGx5KQpgYGAKCgpDb25zaWRlciBhIHNvbWV3aGF0IG5hc3R5LWxvb2tpbmcgY3VydmUgaW4gdHdvIGRpbWVuc2lvbnMuICBXZQpkZWZpbmUgJHokIHRvIGJlIGEgZnVuY3Rpb24gb2YgJHgkIGFuZCAkeSQgdGhhdCBsb29rcyBsaWtlIHRoaXM6IAokJAp6ID0gZih4LCB5KSA9IFxiaWdnbFsxICsgZSBeey18eCArIDJ5fH0gXGNvc1xiaWdnbCgwLjIgXHNxcnR7eF4yICsgeV4yfVxiaWdncikgXGJpZ2dyXSBcdGltZXMgSVx7LTEwIDwgeCx5IDwgMTAgXH0KJCQKSW4gb3RoZXIgd29yZHMsIGl0IGlzIGEgZnVuY3Rpb24gdGhhdCByZXR1cm5zIDAgZXZlcnl3aGVyZSBvdXRzaWRlIG9mIGEgc3F1YXJlCm9mIHdpZHRoIDIwLCBjZW50ZXJlZCBvbiB0aGUgb3JpZ2luLiAgQnV0IHdpdGhpbiB0aGF0IHNxdWFyZSBpdCByZXR1cm5zCnBvc2l0aXZlIHZhbHVlcy4KCkhlcmUgaXMgc29tZSBjb2RlIHRoYXQgZGVmaW5lcyBpdDoKYGBge3J9CnVnbHkgPC0gZnVuY3Rpb24oeCwgeSkgewogIGlmZWxzZSh4IDwgLTEwIHwgeSA+IDEwIHwgeCA+IDEwIHwgeSA8IC0xMCwgCiAgICAgICAgIDAsCiAgICAgICAgIDEgKyBleHAoLWFicyh4ICsgMiAqIHkpIC8gMTApICogY29zKDAuMiAqIHNxcnQoeF4yICogeV4yKSkKICApCn0KYGBgCgpBbmQsIHdlIGNhbiB1c2UgdGhhdCB0byB2aXN1YWxpemUgdGhlIHN1cmZhY2UgYnkgZXZhbHVhdGluZyBhIGxvdCBvZiBwb2ludHMgYW5kCnBsb3R0aW5nIGl0IHdpdGggdGhlIHBsb3RseSBwYWNrYWdlLiAgR28gYWhlYWQgYW5kIHJvdGF0ZSBhbmQgbWFuaXB1bGF0ZSB0aGUgdGhpbmcuCmBgYHtyLCBmaWcud2lkdGg9OH0KdmFscyA8LSBzZXEoLTEyLCAxMiwgYnkgPSAwLjA1KQpwdHMgPC0gZXhwYW5kX2dyaWQoCiAgeCA9IHZhbHMsCiAgeSA9IHZhbHMKKSAlPiUKICBtdXRhdGUoeiA9IHVnbHkoeCwgeSkpCgp6bWF0IDwtIG1hdHJpeChwdHMkeiwgbmNvbCA9IGxlbmd0aCh2YWxzKSwgbnJvdyA9IGxlbmd0aCh2YWxzKSkKCnBsb3RfbHkoeCA9IHZhbHMsIHkgPSB2YWxzLCB6ID0gfnptYXQpICU+JQogIGFkZF9zdXJmYWNlKCkKYGBgCgoKVGhlIHZvbHVtZSB1bmRlciB0aGF0IGN1cnZlIGlzIGdpdmVuIGJ5IHRoZSBkZWZpbml0ZSBpbnRlZ3JhbDoKJCQKXGludF97LTEwfV57MTB9IFxpbnRfey0xMH1eezEwfSBmKHgsIHkpIGR4ZHkgPSBcaW50X3stMTB9XnsxMH0gXGludF97LTEwfV57MTB9IFxiaWdnbFsxICsgZSBeey18eCArIDJ5fH0gXGNvc1xiaWdnbCgwLjIgXHNxcnR7eF4yICsgeV4yfVxiaWdncikgXGJpZ2dyXSBkeGR5CiQkClRoYXQgbG9va3MgYSBsaXR0bGUgaGFpcnkuICBJdCBtYXkgaGF2ZSBhbiBhbmFseXRpY2FsIHNvbHV0aW9uLCBvciBpdCBjb3VsZApiZSBldmFsdWF0ZWQgbnVtZXJpY2FsbHkgd2l0aCBoaWdoIGFjY3VyYWN5LiAgQnV0LCBpZiB5b3UKcmV3cml0ZSB0aGF0IGludGVncmFsIGFzIGFuIGV4cGVjdGlvbiBvdmVyIGEgdHdvLWRpbWVuc2lvbmFsIHVuaWZvcm0gdmFyaWFibGUgb24gdGhlIHNxdWFyZSBmcm9tIC0xMCB0byAxMCBpbiAkeCQgYW5kIGFsc28gZnJvbSAtMTAgdG8gMTAgaW4gJHkkLCB0aGVuIHlvdSBnZXQ6CiQkCjIwIFx0aW1lcyAyMCBcdGltZXMgXGludF97LTEwfV57MTB9IFxpbnRfey0xMH1eezEwfSBmKHgsIHkpIFxmcmFjezF9ezIwfVxmcmFjezF9ezIwfWR4ZHkgPQo0MDAgXHRpbWVzIFxtYXRoYmJ7RX1bZih4LCB5KV0KJCQKd2hlcmUgdGhlIGV4cGVjdGF0aW9uIGlzIHRha2VuIG92ZXIgYSAyLUQgdW5pZm9ybSB2YXJpYWJsZSBmcm9tIC0xMCB0byAxMCBpbiBib3RoIHggYW5kIHkuCgpHaXZlbiB0aGF0IGFzIGEgcHJlYW1ibGUsIHlvdXIgbWlzc2lvbiBpdCB0byBhcHByb3hpbWF0ZSB0aGF0IGludGVncmFsIHVzaW5nIE1vbnRlIENhcmxvIHNhbXBsaW5nLgoKR29vZCBsdWNrISAgVGhlIGFuc3dlciBpcyAkXGFwcHJveCA0MjMuNyQuCgpIZXJlIGlzIGhvdyB0aGF0IGNhbiBiZSBkb25lOgpgYGB7cn0Kc2V0LnNlZWQoMTIzKQp4cyA8LSBydW5pZigxZTYsIC0xMCwgMTApCnlzIDwtIHJ1bmlmKDFlNiwgLTEwLCAxMCkKCmFuc3dlciA8LSAgNDAwICogbWVhbih1Z2x5KHhzLCB5cykpCgphbnN3ZXIKYGBgCg==